// if compiling for Pentium 4 or later uncomment line below
// #define SSE2 
#include "fermiqcd.h"

class open_indices {
public:
  mdp_complex c[4][4][3][3];
};

int main(int argc, char** argv) {
  mdp.open_wormholes(argc, argv);
  define_base_matrices("FERMILAB");

  if(argc<6) {
    mdp << "Usage:\n";
    mdp << argv[0] 
	<< " gaugefile kappa_light csw_light kappa_heavy csw_heavy [ap] [rp] [d1]\n"
	<< "[ap]=20 for 10^-20 absolute precision\n"
	<< "[rp]=12 for 12^-12 relative precision\n"
	<< "[d1]=0.0 rotation amount\n";
    exit(1);
  }

  if(!is_file(argv[1])) {
    mdp << "Unbale to open gauge configuration\n";
    exit(1);
  }
  // read the lattice size from the gauge configuration itself
  mdp_field_file_header header=get_info(argv[1]);
  int* L=header.box;
  int n=3; // number of colors;
  mdp_lattice lattice(4,L,default_partitioning<0>,torus_topology,0,1,false);
  mdp_lattice lattice1d(1,L,default_partitioning<0>,torus_topology,0,1,false);
  gauge_field U(lattice,n);  
  fermi_propagator S_l(lattice,n); // light quark propagator
  fermi_propagator S_h(lattice,n); // light quark propagator
  mdp_site x(lattice); // site on the 4d lattice
  mdp_site x1d(lattice1d); // site on the projected lattice
  int a,b,c,d,i,j,k; // temp variables

  mdp_field<open_indices> open_prop(lattice1d); // field for output open prop

  coefficients light_quark;
  coefficients heavy_quark;

  // making these numbers up
  light_quark["kappa"]=atof(argv[2]);  
  light_quark["c_{SW}"]=atof(argv[3]); 
  heavy_quark["kappa"]=atof(argv[4]); 
  heavy_quark["c_{SW}"]=atof(argv[5]); 
  float ap=1e-20; // inversion absolute precision
  if(argc>6) ap=pow(10.0,-atoi(argv[6]));
  float rp=1e-12; // inversion relative precision
  if(argc>7) rp=pow(10.0,-atoi(argv[7]));
  if(argc>8) light_quark["rotate_field:d1"]=heavy_quark["rotate_field:d1"]
    =pow(10.0,-atoi(argv[7]));
  else light_quark["totate_field:d1"]=heavy_quark["rotate_field:d1"]=0;
  mdp << "Target absolute precision = " << ap << endl;
  mdp << "Target relative precision = " << rp << endl;

  char output_filename[1024];
  sprintf(output_filename,"%s_light_%s_%s_heavy_%s_%s_open",
	  argv[1], argv[2], argv[3], argv[4], argv[5]);

  // load gauge configuration
  U.load(argv[1]);

  compute_em_field(U);
#ifdef SSE2 
  default_fermi_action=FermiCloverActionSSE2::mul_Q;
#endif
  generate(S_l,U,light_quark,ap,rp,5000,rotate_field);
  generate(S_h,U,heavy_quark,ap,rp,5000,rotate_field);

  forallsites(x1d) 
    for(a=0; a<4; a++)
      for(b=0; b<4; b++)
	for(i=0; i<3; i++)
	  for(j=0; j<3; j++)
	    open_prop(x1d).c[a][b][i][j]=0;

  forallsites(x) {
    x1d.set(x(0)); // project a 4d site x into a 1d site x1d
    for(a=0; a<4; a++)
      for(b=0; b<4; b++)
	for(c=0; c<4; c++)
	  for(i=0; i<3; i++)
	    for(j=0; j<3; j++)
	      for(k=0; k<3; k++) {
		d=Gamma5_idx[b];
		open_prop(x1d).c[a][d][i][j]+=
		  Gamma5_val[b]*S_l(x,c,a,k,i)*conj(S_h(x,c,b,k,j));    
	      }
  }
  
  open_prop.save(output_filename);

  mdp.close_wormholes();
  return 0;
}

